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We develop strong-coupling series expansion methods to study two-particle spectra of quantum 
lattice models. At the heart of the method lies the calculation of an effective Hamiltonian in 
the two-particle subspace. We explicitly consider an orthogonality transformation to generate this 
block diagonalization, and find that maintaining orthogonality is crucial for systems where the 
ground state and the two-particle subspace are characterized by identical quantum numbers. We 
discuss the solution of the two-particle Schrodinger equation by using a finite lattice approach 
in coordinate space or by an integral equation in momentum space. These methods allow us to 
precisely determine the low-lying excitation spectra of the models at hand, including all two-particle 
bound/antibound states. Further, we discuss how to generate series expansions for the dispersions 
of the bound/antibound states. These allow us to employ series extrapolation techniques, whereby 
binding energies can be determined even when the expansion parameters are not small. We apply 
the method to the (1+1)D transverse Ising model and the two-leg spin-i Heisenberg ladder. For 
the latter model, we also calculate the coherence lengths and determine the critical properties where 
bound states merge with the two-particle continuum. 
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I. INTRODUCTION 



> 

The study of bound states and multiparticle excitations remains a challenging problem in many-body physics. 
Experimentally, there are several probes for low-dimensional magnetic or strongly correlated electronic systems which 
show spectral features associated with multiparticle continuum and bound states. These include two-magnon Raman 
spectra, optical absorption, photoemission and neutron scattering spectra. The multiparticle features often remain 
poorly understood. On the theoretical side, one example of the intriguing issues that may arise is the role that the 
increasing number of bound states play in the confinement-deconfinement transition in spin-Peierls systems. At the 
transition the spectrum switches from a soliton-antisoliton continuum to elementary triplet excitations, their bound 
states and continuum [jlj. 

A controlled numerical framework for the calculation of multi-particle spectral properties, which can also account for 
various singularities as the coupling constants are varied, is currently missing. In one dimension, a variety of numerical 
, methods including Lanczos, exact diagonalization and most notably density matrix renormalization group (DMRG) 
H hold promise for such calculations. However, unlike ground state and single-particle properties, the calculation of 
. j-h , full dynamical properties like spectral functions still needs more conceptual advances. In higher than one dimension, 
' all of these methods are restricted to small system sizes, which makes it difficult to study the thermodynamic limit. 
On the other hand in the limits of weak or strong couplings, perturbation theory can be used to calculate all 
properties of the multiparticle spectra directly in the thermodynamic limit. If these calculations can be done to high 
orders, one can calculate multiparticle spectra in a systematic manner using extrapolation techniques even in cases 
where the perturbations are not weak. In principle, one can see that as the couplings are increased the number of 
bound states can change and states can come off or merge into the continuum. The resulting singularities should be 
amenable to series expansion methods. 

In this paper, we show how to calculate multi-particle spectral properties from high-order perturbation expansions, 
using a linked cluster method. A brief outline and summary of the work was given in a recent letter H . Our method 
is quite distinct from the flow equation approach of Wegner pi, which has also been used recently by Uhrig and 
collaborators |^,^] for the study of multiparticle spectral properties in one and two dimensions. 
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The linked cluster method is by far the most efficient way to generate perturbation series expansions for quantum 
Hamiltonian lattice models. For the ground state energy and related properties, a linked cluster approach was first 
discussed in unpublished work by Nickel M, followed by work of Marland M, Irving and Hamer M and others, as 



reviewed by He et al 10 1. The approach was later rediscovered and applied to a whole new range of problems in 



condensed matter physics by Singh, Huse and Gelfand |ll| , |l2j |. 

For the energies of excited states, it is more difficult to formulate a true linked cluster expansion, although related 
methods have been known for some time |Q,|l3|- It was only in 1996 that the key to a true linked cluster expansion 
for one-particle excited states was discovered by Gelfand [ fhij . Since then, many applications of this technique have 
been made, calculating single-particle energies, dispersion relations and spectral functions in models of interest in 
condensed matter physics. For a recent review, see Gelfand and Singh p"5| . 

The key advantage of the cluster expansion method is that the calculations can be carried out systematically and 
efficiently by fully automated computer programs. Furthermore, these methods work by breaking up the thermody- 
namic problem into a purely combinatorial problem and a number of finite-cluster problems. Thus, while they are 
technically harder in higher than one dimension, the difficulty is not fundamental. In fact, over the years, a num- 
ber of workers |jl^ , ^ , ^5|Jl^Jl9[] have independently developed efficient computer programs to generate these clusters 
automatically, and the cluster data up to quite high number of vertices for most 2-dimensional and 3-dimensional 
lattices including the simple-cubic, BCC and FCC lattices have been generated fl6|| : these data can be applied to a 
wide range of models. 

At the heart of our new approach is a generalization of Gclfand's linked cluster expansion for single-particle excited 
states to two-particle states. From a technical point of view, our most notable achievement is the development of 
an orthogonality transformation which leads to a linked cluster theorem for multi-particle states even when their 
quantum numbers are identical with the ground state. We show how to calculate energies and dispersion relations for 
two-particle excitations, and coherence lengths for the bound states. The further generalization to higher number of 
particle is then obvious in principle. 

As a first check to ensure that the method is working correctly, we apply it to the case of the transverse Ising model 
in one dimension, which can be solved exactly in term of free fermions [|o| . We show that the series for the 2-particle 
state agree with the exact results up to 12th order. 

Finally, we apply the method to a non-trivial model, the 2-leg spin-i Heisenberg ladder, which has been much 
discussed recently pl|-p8| as a prime example of a one-dimensional antiferromagnetic system with a gapped excitation 
spectrum. The two-particle bound states have already been studied by Damle and Sachdev [^5| and Sushkov and 
Kotov [ p7p^ ]. We perform a detailed study of these bound states, exhibiting in particular the characteristic features 
as each bound state emerges from the continuum. 

In a companion paper p9[ | , we apply the same techniques to a still more interesting case, the frustrated alternating 
Heisenberg chain, which displays the confinement-deconfinement transition discussed by Affleck jlj . 

The organization of the paper is as follows. Section II of the paper lays out the formalism and methods used to 
obtain the 2-particle spectra. Section III discusses the applications to the transverse Ising chain and the Heisenberg 
ladder, and Section IV summarizes our conclusions. 



II. FORMALISM 

We consider a Hamiltonian 

H = H + \H ll (1) 

where the unperturbed Hamiltonian Hq is exactly solvable, and A is the perturbation parameter. In the lattice models 
of interest here, Hq will typically consist of single-site operators, while interaction terms between different sites will 
be included in the perturbation operator Hi. The aim is to calculate perturbation series in A for the eigenvalues of 
H and other quantities of interest. The calculation proceeds in three stages. 



A. Block Diagonalization 

On any finite lattice or cluster of sites, the first step is to 'block diagonalize' the Hamiltonian to form an "effective 
Hamiltonian" , where the ground state sits in a block by itself, the 1-particle states form another block, the 2-particle 
states another block, and so on. Here a "particle" may refer to a lattice fermion, a spin-flip, or other excitation, 
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depending on the model at hand. We assume that all the unperturbed states in each block are degenerate under Ho. 
There is no unique way to block diagonalize the Hamiltonian, but the eigenvalues and final results should be unique, 
independent of the method used, as long as the cluster expansion works correctly. Gelfand E3] used a similarity 
transform for this purpose: 

H eS = q-1 H q ( 2 ) 

This works correctly for most 1-particle problems, and also for those 2-particle states which have different quantum 
numbers to the ground state. However in general, especially for the 2-particle states which have identical quantum 
numbers to the ground state, we need to be a little more careful than this, in order to preserve all the proper 
symmetries of the Hamiltonian. We must ensure that the transformation is unitary. Here we will only consider the 
case when the Hamiltonian is real symmetric, and can be block diagonalized by an orthogonal transformation, 

H cS = T HO (3) 

or more conveniently 

OH cS = HO , (4) 

where 

O t = O- 1 . (5) 

The orthogonality of O can be ensured by writing 

O = e s , (6) 

where S is real, antisymmetric 

S T = -S . (7) 

This transformation is constructed order-by-order in perturbation theory. The matrix elements of H between 
different blocks are zero, up to the given order in perturbation theory. Each matrix is expanded in powers of A: 

oc 

= X x "° n) > ( 8 ) 

n=0 



S = xnsn) > ( 9 ) 



n=0 



# cff = E A "Cff> (io) 



n=0 



where at zeroth order we set 



S°)=Q, 0°)=/, H°J S = H Q . (11) 



where I is an unit matrix, Hg is a diagonal matrix, with diagonal matrix elements Ef . 
At higher orders n ^ 0, we have 



n 1 n 

0«) = S™) + - J2 S m) S l) 6 m+lin + - £ S m ^S l '>S k ^ m+l+k ^ + ... (12) 

m.l—1 m,l,k—l 

and 

n 

X O m ^H%5 m+Ln = H O n ^ + i/xO"- 1 ) (13) 



m,l=0 
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FIG. 1. Block structure of the matrices H™g and S n ' . Setting the upper right blocks of H™g to zero determines the corre- 
sponding (shaded) blocks of S 1 "'; the diagonal blocks of S"^ are set to zero. 



and it is convenient to define 

R n) = Q n) _ gn) 

If we demand that at any given order n the off-diagonal blocks of H cS in (say) the upper right triangle vanish, 
then equations (13) determine the entries in the corresponding blocks of (Figure [l]). The transposed blocks in the 



lower left triangle are then determined by the antisymmetry condition (|?|) ; and only the diagonal blocks of 5* remain 
to be determined. The simplest choice is to set the diagonal blocks to zero. Thus S 1 ™-* is completely determined: 

n— 1 

S S = - R $ + (E o_ EO) { H i° n ~ 1) E O^H%5 m+Un } l0 (15) 

3 1 ' mj=l 

or 

71—1 

°t = JW^W) {Hl0n ^ E O^H^S^nhj (16) 

^ 3 i' m,l=l 

for elements ij in the off-diagonal (shaded) blocks. Then 



n-l 



(Ksh = - O m) ^«m+i,n}« (17) 

m,l— 1 



for elements in the diagonal blocks. The right-hand sides of equations ( |15| , |16|]17| ) can all be computed from the results 
at order (n— 1). 

The key differences here from the similarity transformation are as follows. In the similarity transformation, the 
diagonal blocks of O n ' are undetermined, and so are chosen to be zero, while the off-diagonal blocks of O n ^ are 

n) 

antisymmetric and can be determined by demanding the off-diagonal blocks of H c ^ to be zero. In the orthogonal 
transformation, on the other hand, the diagonal blocks of O™' cannot be chosen to be zero. Instead the diagonal 
blocks of S n ' are chosen to be zero, while the diagonal blocks of O n ^ are required to be nonzero by orthogonality, and 
are determined by Eq. ([TJ]). 

At the end of this process, the effective Hamiltonian has been block diagonalized, up to a given order in perturbation 
theory. The orthogonal transformation will transform the unperturbed two-particle states into "dressed" states 
containing admixtures of different particle numbers; and in particular, there will be no annihilation process for these 
"dressed" states. The states will still be labelled by the positions of the original unperturbed particles; but now they 
will contain admixtures of other particle states at nearby locations. 

At any finite order in perturbation theory, we may assume that the effective Hamiltonian will remain "local" (that 
is, interactions between states will not extend beyond a finite range); and will have the same bulk symmetries as the 
original Hamiltonian, such as translation symmetry. These properties are sufficient to admit a linked cluster approach 
to the calculation of eigenvalues. 

We note that the solution of the equations above is not nearly as efficient as the similarity transformation of Gelfand: 
in particular, the solution of the equation ( |T^ ) is expensive in CPU time and memory. In the Appendix, we discuss an 
alternative '2-block' scheme which has the same efficiency as Gelfand's; but which does not always allow a successful 
cluster expansion. 



B. Linked Cluster Expansions 

Let us briefly summarize the linked cluster approach in various sectors. 
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FIG. 2. Decomposition of a 1-particle matrix element into irreducible components. The round box denotes the full matrix 
element, the square boxes the irreducible matrix elements, and the single line denotes a delta function. 



1. Ground-state energy 



The ground-state energy Eo is a simple extensive quantity, and obeys the "cluster addition property" : if C is a 
cluster (or set of sites and bonds on the lattice) which is composed of two disconnected sub-clusters A and B, then 



EC=e£ + E*. (18) 

Hence one finds @-|l2| that the ground-state energy per site for the bulk lattice can be expressed purely in terms of 
contributions from connected sub-clusters a: 

eo = y^J a e a , (19) 

a 

where l a is the "lattice constant" , or number of ways per site that cluster a can be embedded in the bulk lattice, and 
e a is the "proper energy" or "cumulant energy" for the cluster a. In the language of Feynman diagrams, e a can be 
thought of as the sum of all connected diagrams spanning the cluster a [ |30| |. 

A similar formula holds for the ground-state energy of any connected cluster a with open boundaries: 

E^ = J2C^ , (20) 



where Cp is the embedding constant of the connected sub-cluster (3 within cluster a. 

Equations (19) and ( |20| ) form the basis for a simple and efficient recursive algorithm to generate a perturbation 
series for e - The steps are: 

i) Generate a list of clusters a, with their lattice constants l a and embedding constants C^, appropriate to the 
problem at hand jiJ0|l7]]; 

ii) For each cluster a, the diagonal entry in the 0-particle sector of H cS gives a perturbation series for the energy 

iii) Now invert equations (|2(]) to solve for the cumulant energies e Q , and substitute in (|l9|) to obtain the desired 
perturbation series for eo- 



2. 1-particle excited states 

Gelfand Jl4|] discovered how to generalize the approach above to one-particle excited states. Let 

i? 1 (i,j) = (j|^ off |i) (21) 

be the matrix element of H cS between initial 1-particle state |i) and final 1-particle state labelled according to their 
positions on the lattice. The excited state energy is not extensive, and does not obey the cluster addition property; 
but there is a related quantity which does. If cluster C is made up of disconnected sub-clusters A and B, and states 
|i) and |j) reside (say) on cluster A, then 

J B 1 c (i,j)=^(i,j)+ii; B . (22) 
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But if we define the "irreducible" 1-particle matrix element (Fig. ||) 

Ai(i,j)=£?i(iJ)-£b*j , 

then 



Af(i,j) = Af(i,j), 



whereas if |i) and |j) reside on cluster B, then 



or in general 



Af(i,j) = Af(iJ) 



Af(i,j) = Af(i,j)+Af(i,j) 



(23) 



(24) 



(25) 



(26) 



where Ai(i,j) vanishes for any cluster not containing i and j. Note that a 1-particle state cannot annihilate from one 
sub-cluster and reappear on the other, after the initial block diagonalization. 

From the cluster addition property ( p6[ ) it follows that the elements Ai(i, j) can be expanded in terms of contribu- 
tions from connected clusters alone, which are also "rooted", or conne cted to the positions i and j. Hence they can 
be calculated efficiently by an algorithm like that of subsection ( [IB 1 ). 



3. 2-particle states 



The generalization to two-particle states is now not hard to find. Let 

effi 



£ 2 (i,j;k,l) = (k,l|ff on |i,j 



(27) 



be the matrix element between initial 2-particle state |i,j) and final state |k, 1). To obtain a quantity obeying 
the cluster addition property, we must subtract the ground-state energy and 1-particle contributions, to form the 
irreducible 2-particle matrix element [Fig. || : 



A 2 (i, j; k, 1) = E 2 (i, j; k, 1) - E (5 iM 6 iA + S M <5j, k ) - Ai(i, k)<5 ja - Ai(i, l)dj,, k 
-Aiakjay-Aia,!)**. 



(28) 



This quantity is easily found to be zero for any cluster unless i, j, k and 1 are all included in that cluster, and it obeys 
the cluster addition property. Once again, the block diagonalization ensures that two particles cannot "annihilate" 
from one cluster and "reappear" on another disconnected one. Thus the matrix elements of A 2 can be expanded in 
terms of connected clusters alone, which are rooted or connected to all four positions i, j, k, 1. 






FIG. 3. Decomposition of two identical particle matrix element into irreducible components. Notation as in Fig. 2. 



C. Calculation of eigenvalues 



For the ground state energy, a perturbation series for the eigenvalue was already obtained at the end of subsection 
[IB 1. For the excited state sectors, some further work is required. 
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1. 1 -particle states 

A perturbation series for the dispersion relation of the 1-particle states can be calculated by a Fourier transform. 
Translation invariance implies that 

A 1 (i,j) = A 1 (6), (29) 
where 5 is the difference between positions i and j; and that the 1-particle states are eigenstates of momentum: 

i K > = -^E ex p( iK -j)ij) > ( 3 °) 
j 

(where N is the number of sites in the lattice), with energy gap 

wi (K) = ^ Ai {5) cos(K • 6) . (31) 

<5 

Here we have assumed that Ai(J) is reflection symmetric, so that 

A 1 (-5) = Ax(<5) . (32) 

2. 2-particle states 

The calculation of the eigenvalues in this case is a little more involved than in the 1-particlc case. We follow the 
procedure of Mattis Q . 

Consider an unsymmetrized state of non-identical particles, types a and b. Then there are N(N — 1) states on an 
TV-site lattice, labelled by positions |i,j), where i,j refer to the positions of particles a and b, respectively. We have 
assumed here that two particles may not reside at the same position (the results are easily amended if this is not the 
case). Then the irreducible 2-particle matrix element is: 

Af(iJ;k,l) = £f(i,j;k,l) 

-£ <5i,iAi - A?(i, k)*j,i - A$(j, I)5 ilk . (33) 



Let the 2-particle eigenstate be: 



substitute in the Schrodinger equation 



E/yliJ), (34) 



H cB \^ =E\i>) (35) 
and take the overlap with (i,j|, then one obtains: 

(S-^/y-^A^i^kj-^A^kJ^^^AfCk,!;^^!, (i^j). (36) 
k#j k#i k,i 

Completing the sums on the left-hand side, one obtains: 

(E - E )f u - X)[A?(k, i)/kj + A?(k, j)/ ik ] = J2 A^ b (k, 1; i, j)/ kl 

k k,l 

-A;a,i)/jj-Aj(i,j)/ u , (37) 
The fictitious amplitudes fa are introduced to simplify the calculations, and are taken to be defined by these equations 
Now define a centre-of-mass position co-ordinate 
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and relative co-ordinates 



R= J(i+j + k + i) , 



r = i(i + j-k-I) 



<*i = i-j , 

<y a = k - 1 . 

Translation invariance then implies that 

A 2 (iJ;k,l) = A 2 (r,S 1 ,5 2 ) . 

Next, perform a Fourier transformation, 

^ K ' f i) = ^E el(kl ' i+k2 ' j) ^' 

u 

where K, q are the centre-of-mass and relative momenta: 

K = (ki +k 2 ) , 

q = \( k i -k 2 ) , 

then equation (|37]) leads to 

[E-Eq- 5^t A ?(*) cos ( K ' <5/2 + q • 5) + cos(K ■ 5/2 - q ■ <5)]/(K, q) 

^^/(K,q')[ J] Af(r,<5 1 ,<5 2 )cos(KT + q. ( 5 1 -q'.5 2 ) 

- ^ A?(<$) cos(K • 5/2 + q • 5) + Aj (<S) cos(K • <5/2 - q • 5)] , 

where we have again assumed reflection symmetry 

Ag 6 (r,fr,&) = A3 6 (-r, -*!,-&). 
Finally, look for solutions with definite exchange symmetry. 
Symmetric states 

therefore 

/(K,q) = +/(K,-q). 
'Averaging' over /(K,±q) (i.e. taking i[/(K,q) + /(K,-q)]), we get: 

{£ - B - Et A ?W + A iW] C0S ( K ' <V 2 ) cos (q ' <W( K ' ^ = 



l^/(K,q')[ ^ Af(r,<5 1 ,<5 2 )cos(K.r)cos(q- ( 5 1 )cos(q'- ( 5 2 ) 
- 5^[A?(<y) + Aj(5)] cos(K • 5/2) cos(q • (5)] . 



AT 

q' r,oi,02 



Antisymmetric states 



fij — fji (52) 

therefore 

/(K,q) = -/(K,-q). (53) 
'Averaging' over /(K, ±q) (i.e. taking |[/(K,q) — /(K,— q)]), we get: 

{E-E - £[A?(*) + A*(<J)] cos(K • ,5/2) cos(q • <5)}/(K, q) = 

s 

l^/(K,q') £ Af(r^ 1 ^ 2 )co S (K-r)sin(q- ( 5 1 )sm(q'^ 2 ). (54) 

Identical particles 

If the particles a and b are identical, the solution is the same as for symmetric states except the labels a and b must 
now be dropped, and to avoid double counting it turns out that the A2 term must be multiplied by an extra factor 
of 1/2: 

[E - E - 2 A i( S ) C0S ( K ' <V 2 ) cos(q • 5)}f(K, q) = 

s 

l^/(K,q')[i A 2 (r,«5 1 ,(5 2 )cos(K.r)cos(q. ( 5 1 )cos(q'. ( 5 2 ) 

q' r, 6 1,5-2 

-2 Ai(<5) cos(K • 5/2) cos(q ■ 6)] . (55) 
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The above integral equations can be solved, for a given value of K, using standard discretization techniques. Instead 
of using continous momentum q, one can use N discretized and equally spaced values of momentum, so that instead 
of solving the complicated integral equation, one only needs to compute the eigenvalue and eigenvector of an N x N 
matrix for the discretized system. Notice that the matrix is nonsymmetric due to the unphysical fa term we have 
introduced in Eq. (|37]), but even so the eigenvalues obtained from this matrix are real. The solutions we obtain 
also include an unphysical one with eigenvalue equal to (this is also due to the unphysical fa term). The results 
obtained from the calculation with discretized momenta will converge to those with continous momentum as N — > 00. 
Actually for those bound states with finite coherence length, the calculation will normally be well converged for quite 
small values of N, but for unbound states, we have an infinite coherence length, so one may need to do finite N 
extrapolations to get results at N = 00. 

There are two methods to compute the eigenvalues of the matrix for the discretized system. Obviously one can 
get numerical results for the eigenvalues, for a given value of coupling A and momentum K, via standard numerical 
techniques where we just perform a naive sum for the series in Aj and A 2 . The results presented in a preceding Letter 
H are based on this method; but then one cannot carry out a series extrapolation, and so one may not be able to reach 
a region of critical coupling. A better technique is to compute the series in A for the eigenvalues through degenerate 
perturbation theory: that is, by explicit diagonalization of the matrix within the degenerate subspace, order by order 
in perturbation theory, and then one can perform a series extrapolation. The problem with this method is that the 
series does not always exist, for example for those bound states appearing at some nonzero value of A. 

The two particle continuum is delimited by the maximum (minimum) energy of two single particle excitations 
whose combined momentum is the center of mass momentum. Apart from the unphysical eigenvalue, there may be 
multiple solutions above/below the two-particle continuum. Those solutions with energy below the bottom edge of 
the continuum are the bound states, while the solutions with energy higher than the upper edge of the continuum are 
the antibound states. The binding energy is defined as the energy difference between the lower edge of the continuum 
and the energy of the bound state, while the antibinding energy is defined as the energy difference between the upper 
edge of continuum and the energy of an antibound state. 

Note that the series for A 2 may depend on the transformation used to block diagonalize the Hamiltonian. If we 
compute A 2 (and also Ai) to order n, the resulting series for the 2-particle energy obtained from the above integral 
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equation will have two parts: the part up to order n is independent of the transformation, while the higher order 
terms are incomplete, and may depend on the transformation. The numerical solution of the integral equation may 
also depend partly on the transformation, since it contains the higher order term. Also note that the series for A2 
need not have any singularities. The singularities, if they exist, arise in the solution of the Schrodinger equation, so 
our method should be able to explore new bound states arising as we vary the momentum K. If we get a numerical 
solution, rather than a series solution, to the Schrodinger equation, we should also be able to explore new bound 
states arising as we increase A as long as the naive sum to the series converges. 

D. Finite Lattice Approach 

Once the cluster expansions for the irreducible matrix elements Ai and A2 have been developed, the Schrodinger 
equation in the two-particle subspace can be solved by an alternative method that works in coordinate space rather 
than momentum space. By restricting to a finite but large system with periodic boundary conditions, the two-particle 
Schrodinger equation becomes a finite-dimensional matrix of equations. The cluster expansion results provide the 
matrix elements of the effective Hamiltonian as a power series in the expansion parameter. The centre of mass 
momentum is a conserved quantity, thus, for a given value of the centre of mass momentum, one is left with a 
Schrodinger equation in the separation variable. One can truncate the perturbation theory at a given order and solve 
the Schrodinger equation numerically. One can then vary the size of the system, which only increases the dimension 
of matrix to be diagonalized linearly, to study convergence. We have frequently used this method to compare with 
and check the momentum-space discretization solutions. 

This 'finite lattice approach' also allows us to obtain power series expansions for bound state energies, by a non- 
degenerate perturbation theory, provided the bound state exists "localized" in the limit A — > 0. For those "extended" 
bound states in the limit A — > one still needs to do degenerate perturbation theory, just as in the case of 
the momentum-space discretization solutions. Another advantage of this method over the the momentum-space 
discretization technique is that the matrix one deals with is always symmetric. 

Furthermore, it gives us explicit real-space wave functions, from which the coherence length and other properties 
can be deduced. The coherence length L is defined by 

V Irlf 2 

L = Ml (56) 
where / r is the amplitude (the eigenvector) for two single-particle excitations separated by distance r. 

III. RESULTS 

We apply the new method to the (1+1)D transverse Ising model and a two-leg spin-i Heisenberg ladder. 

A. Transverse Ising model 

In order to verify that our new technique is giving the correct results, we firstly apply it to a simple model, the S = \ 
transverse Ising model in (l+l)-dimensions, which is exactly solvable in terms of free fermions. The Hamiltonian for 
it reads 

# = £(!-<) -A^?^, (57) 

i i 

Here we take the first term as the unperturbed Hamiltonian H$, and the second term as the perturbation H\. The 
ground state of Ho is the unique state with all spins pointing up. The lowest excited states (1-particle excitations) 
for Ho flip one of the spin from spin up to spin down. The exact result |fjO| for the 1-particle dispersion relation is 

Ei(q) = 2Vl + A 2 -2Acos<7. (58) 

For the 2-particle excitations, the unperturbed states have two spins down. Since this model can be mapped into 
free fermions, there are no 2-particle bound states, and the 2-particle excitation energy is simply the sum of two 
1-particle dispersions, that is 



10 




FIG. 4. The binding energy (full points) and antibinding energy (open points) Et versus 1/N 2 (N is the size of the matrix) 
for the transverse Ising model with coupling A = 0.5 and momentum k — (dotted lines), 7r/2 (dashed lines), n (solid lines). 



E 2 {qi, <h) = 2\/l + A 2 - 2Acosgi + 2y/\ + A 2 - 2Acosq 2 , (59) 

where q\ and qi are the momenta of each particle. Note that this is a non-trivial example for our method as the 
similarity transformation does not even lead to a cluster expansion. 

We have implemented the algorithm described above for this model. For the 1-particle excitation, we can easily 
reproduce the exact results through the different block diagonalization schemes mentioned before. For the 2-particle 
excitations, although there are no bound states, the terms A 2 are not zero. That is because we are using the spin 
representation; in a fermion representation, the A 2 would be expected to vanish. We have computed them to order 
A 12 by using the 2-block method. The series coefficients up to order A 6 are given in Table Q. With these series, 
one can solve the discretized version of the integral equation to get the binding and antibinding energy for any given 
value of momentum k and coupling A. Our results show that for all k and A, the binding/ antibinding energy scales 
as 1/N 2 , and approaches to zero as N — > oo: this is consistent with the absence of bound/antibound states in this 
model. The results for A = 0.5 and k = 0, 7r/2, tt are shown in Fig. |4|. We have also checked that the resulting series 
for E2 agrees with ([59]) for the lowest and highest energy of 2-particle states up to 12th order, and the coherence 
length is infinity, as expected. 



B. Heisenberg ladder 

The second model we have investigated is the 2- leg spin-i Heisenberg ladder, where the Hamiltonian is 

H = ]T{ J±S, ■ S * + J t S < ' S *+i + S * ' S ' i+ i]} > (60) 

i 

where Sj (SQ denotes the spin at site i of the first (second) chain. J is the interaction between nearest- neighbor 
spins along the chain, and Jj_ is the interaction between nearest-neighbor spins along the rungs. In the present paper 
the intra-chain coupling is taken to be antiferromagnetic (that is, J± > ) whereas the interchain coupling J can be 
either antiferromagnetic or ferromagnetic. 

The antiferromagnetic Heisenberg ladder has attracted a good deal of attention recently |2l|-|2^| . It is of experimental 
interest in that there are a number of quasi-one-dimensional compounds which may be described by the model |pl| . 
It is also a prime example of a one-dimensional antiferromagnetic system with a gapped excitation spectrum. Damle 
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FIG. 5. The excitation spectrum for the Heisenberg spin ladder at J/J± = 0.2. Beside the two-particle continuum (gray 
shaded), there are three massive quasiparticles: a singlet bound state (S), a triplet bound state (T) below the continuum and 
a quintet antibound state (Q) above the continuum. 



and Sachdev |2J| as well as Sushkov and Kotov |27]] have shown that the system exhibits two-particle bound states, 
one singlet and one triplet. Our aim here is to explore the properties of these bound states more closely. 

In the dimer limit J/Jj_ = 0, the ground state is the product state with the spins on each rung forming a spin 
singlet. The first excited state consists of a spin triplet excitation on one of the rungs. As J/J± increases, this 
state evolves smoothly, and the system has a gapped excitation spectrum The dimer expansions have been 

computed previously up to order (J/J±) 23 for the ground-state energy and up to order (J/J±) 13 for the 1-particle 
triplet excitation spectrum j23) . The occurrence of two-particle bound states in this model has been shown by first- 
order strong-coupling expansions pq,E6| as well as a leading order calculation using the analytic Brueckner approach 



Here we have calculated series for the dispersions of the 2-particle bound states up to order (J/J±) 7 for the singlet 
bound state (S), and to order (J/J±) 12 for the triplet bound state (T) and the quintet antibound state (Q). The 
reason why the singlet series is computed to only 7th order compared to 12th order for the triplet and quintet states is 
that the singlet has the same quantum numbers as the ground state. Thus a much more elaborate orthogonalization 
method is required to implement the cluster expansion for the singlet. For the triplet and quintet bound states, we can 
use the similarity transformation or the 2-block orthogonal transformation to implement the cluster expansion. Up to 
order (J/J±) 3 , the dispersion for the singlet bound state (Es/J±), triplet bound state (Et/J±), quintet antibound 
state (Eq/J±) are 

3x 19 x 2 9x 3 fx x 2 51x 3 \ 

E s /J± = 2 1 1 h 1 cos(fc) 

b/ ^ 2 16 32 V 2 8 128 J y ' 

5x 2 21x 3 \ , , N 37 a; 3 cos(3 k) „,, s . . 

~~L6 32~) C0S(2 k) 128^ + °^ < (61&) 

„ , T „ 3x 11 x 2 17 x 3 ( x 2 9x 3 \ 
*Wx=2- T + — + — +(-*- T + — j cos(fc) 

(61b) 

3x llx 2 3x 3 ( x 2 27x 3 \ 
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FIG. 6. The binding/antibinding energy Et/J± (upper window) and the rescaled bindhig/antibinding energy Et/ J±(k — k c ) 2 
(lower window) for the Heisenberg spin ladder at J/J± = 0.2. 



/ x 2 3x 3 \ 7a; 3 cos(3fc) _ 4s , . 

+ 1-T-T) cos(2fc) + 16^ + ° ( } ' (6 c) 

where x = J/Jj_. The full dispersion series for the singlet state, and the series for the energy gap at k — tt for singlet, 

tlet and quintet bound/ antibound states and the lower edge and upper edge of continuum are listed in Table [n| and 
the other series are available upon request fl3^] . Figs. |f^and |g] show the dispersion and the binding/antibinding 
energy at J/Ji = 0.2 for the two-particle continuum as well as the the two-particle bound/antibound states. Here we 
can see there is a singlet (S = 0) and a triplet (S = 1) bound state of two elementary triplets below the two-particle 
continuum, and a quintet (S = 2) antibound state above the continuum. 

From these graphs, we can also see that these bound / antibound states exist only when the momentum k is larger 
than some "critical momentum" k c : the series in Eq. ( plj) and Table |l| are valid only for k > k c . It is interesting to 
explore the behaviour of the binding energies near this critical momentum. From the series for the one-particle and 
two-particle dispersions, one can get leading order results for k c as 

C VW^+0(x) _ 5 = 

k c = I 2tt/3 - 5x/(2V3) - 109x7(48^/3) + 0{x 3 ) S = 1 (62) 
{ 2tt/3 + 5x/(2\/3) + 47x 2 /(48V3) + 0(x 3 ) S = 2 

and in the limit k — > k c , the behaviour of the binding energy is 

E b /J = (k — k c ) 2 [5x/8 + 975a; 2 /128 + 0(x 3 )] 

+(k - k c ) 3 [12 + 115a: + 0(x 2 )]Vmi/192 + 0[(k - k c f] (63) 

for the singlet bound state, and 

E b /J =(k- fc c ) 2 [3/8 - x/32 + 0.45313x 2 + 0(x 3 )} 

+ {k - /c c ) 3 [\/3/16 - 53x/(64\/3) + 0.19245a; 2 + 0{x 3 )} + 0[(k - fc c ) 4 ] (64) 

for the triplet bound state. For the quintet antibound state, the antibinding energy is 

E b /J ={k- k c ) 2 [3/8 + x/32 - 0.40625x 2 + 0{x 3 )] 

+ (k - fc c ) 3 [V3/16 + 53x/(MV3) + 1.00886x 2 + 0{x 3 )] + 0[(k - k c f] . (65) 
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FIG. 7. The critical (J/J±) c versus k for singlet (S), triplet (T) and quintet (Q) bound/antibound states of the Heisenberg 
ladder. The solid lines with errorbars are the results obtained from Dlog Pade approximants to the series for binding/antibinding 
energy Eb for given k, while the dotted lines are the results of Eq. (p2[). 



Here one can see that for all bound/antibound states the "critical index" is 2, independent of the order of expansion, 
so one expects that this is exact. 

A better way to locate the critical line in the (J/J±)-k plane is to calculate the Dlog Pade approximants to the 
series for the binding/antibinding energy at a fixed momentum k. For those critical points lying at x c < 0.2, the 
resulting critical point and critical index are very accurate, correct up to 5 digits, and again one finds the critical index 
is exactly 2. The results for the triplet bound state at k = Sir/ 5 are given in Table IV. The results for the critical 
points are given in Fig. |t], together with the results from Eq. ([32]). From this figure, one can see that as J/Jx — * °o, 
k c for the singlet and triplet bound states approaches the same value, about 0.47T, while k c for the quintet antibound 
state approaches it. To demonstrate that Eb is proportional to (k — k c ) 2 near k c , we also plot in Fig. [| the results for 
E b /J±(k - k c ) 2 at J/J ± = 0.2. 

The binding/antibinding energy at k = tt for bound/antibound states versus J/J±_ is plotted in Fig. [|. In the 
limit J/J± — ► 0, Eb/J± is proportional to J/J±, so in the figure we plot Et/J versus J/J±. We can see that as 
J/J± increases, E b j J for the singlet bound state firstly increases, passes through a maximum at about J/J± = 0.4, 
then decreases, while E^/J for the triplet bound state and the quintet antibound state decreases monotonically. At 
J I Jj_ = 1/2, we find the binding/antibinding energies at k = it for the singlet, triplet, and quintet bound/antibound 
states are Et/J =1.03(3), 0.385(1) and 0.0855(5), respectively. The binding energy for the singlet bound state is 
substantially larger than the value 0.70 obtained in [p8| . 

We have also computed the coherence length L for these bound/antibound states. The results for J/J± = 0.2 are 
shown in Fig. where we find that L diverges as l/(fc — k c ) as k approaches k c . This is to be expected, as the state 
becomes unbound at that point. The coherence length at k = it versus J/J± is shown in Fig. [lO] where we can see 
that at J = 0, L = 1. This is as expected, as the formation of these bound states is due to the attraction of two 
triplets on neighboring sites. As J/J± increases, the coherence length L increases slowly. L for the quintet antibound 
state is larger than that for the triplet bound state, which is larger than for the singlet bound state. 



IV. CONCLUSIONS 



In conclusion, we have developed new strong-coupling expansion methods to study two-particle spectra of quantum 
lattice models. We described in full detail the block diagonalization of the Hamiltonian, order by order in perturbation 
theory, to construct an effective Hamiltonian in the two-particle subspace. This work is closely related to Gelfand's 



14 




0.2 0.4 0.6 0.8 1 



J/J, 

FIG. 8. The binding/antibinding energy Eh at k — it versus J/Jx for singlet (S), triplet (T) and quintet (Q) bound/antibound 
states of the Heisenberg ladder. Several different integrated differential approximants to the series are shown. 
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FIG. 9. The coherence length L versus momentum k for singlet (S), triplet (T) and quintet (Q) bound/antibound states of 
the Heisenberg ladder at J/J± = 0.2. 
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FIG. 10. The coherence length L versus J/J± at k — it for singlet (S), triplet (T) and quintet (Q) bound/antibound states 
of the Heisenberg ladder. 



prior work on using similarity transformations to obtain an effective Hamiltonian for the single-particle subspace 
. We found that one needs to define a two-particle irreducible matrix element, for a cluster expansion to exist. 
Furthermore, one needs to maintain explicit orthogonality in the transformations in order to study the two-particle 
subspace characterized by identical quantum numbers to the ground state. An example of the latter is the two-particle 
singlet excitation sector in dimerized spin models. 

We have discussed the solution of the integral equation one obtains by a Fourier transformation of the two-particle 
Schrodinger equation and by a 'finite-lattice approach'. These allow us to precisely determine the low-lying excitation 
spectra of the models at hand, including all two-particle bound/antibound states. Furthermore, we have shown that 
one can generate series expansions directly for the dispersions of the bound/antibound states, provided these bound 
states exist in the limit A — > 0. These allow us to apply series extrapolation techniques such as Dlog Pades and 
differential approximants to study binding energies even when the perturbation parameter is not small. 

We applied the method to the (1+1)D transverse Ising model and the two-leg spin-i Heisenberg ladder. While 
the first model does not include any bound states, we find a singlet and a triplet bound state in the latter model as 
well as a quintet antibound state. We generated explicit expressions for the dispersions of these states as series in 
the exchange couplings. Further, we have determined the critical momenta k c , where these additional massive quasi- 
particles merge with the two-particle continuum, which are non-zero for all three states. The explicit expressions of 
the binding energies at the respective critical momenta are found to contribute first in order (k — k c ) 2 , independent 
of the order of the strong coupling expansion. We computed the coherence length for these states and find that the 
coherence length diverges as one approaches the critical momentum where these states become unbound. 

There are several possible direction for future research along these lines. Of course there are many different models 
to which these methods might be applied. In particular, it remains to show that the linked cluster expansion works 
sucessfully for two or higher dimensional models. One would also like to know how to calculate other quantities 
associated with multipartcle excitations, such as spectral weights, lifetimes, and scattering S-matrices. The latter 
would provide a handle on some important dynamical properties of the system. 
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APPENDIX A: 2-BLOCK APPROACH 



There is an alternative way to perform the block diagonalization of Section II A , which is almost as efficient as 
Gelfand's similarity transformation. The idea is to separate the effective Hamiltonian into only two blocks, one 
containing the states in the sector of interest (e.g. the 1-particle states, or the 2-particle states), and the other 
containing all other states. One can prove in this two-block approach that O"' determined in this way is antisymmetric 
with respect to the off-diagonal blocks, and symmetric with respect to the diagonal blocks. Rather than use the 
complicated equation ( |12|) , one can then determine the diagonal blocks of O n ' in a much more efficient way by the 
orthogonality condition (pO which can be rewritten in the following form: 

n-l 

{O n) + O n)T } l3 = - {O m) O n - m)T } l3 (Al) 

m=l 

for elements in the diagonal blocks. Thus one can dispense with the matrix S, and work with O only. 

Unfortunately, although it is more efficient, this approach does not always seem to allow a successful cluster 
expansion. The reason for this is not understood at the present time. 
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TABLE I. Series coefficients for A2(r, Si, 82) = ^ fe A§ (r, 81, S 2 )x k in the (1+1)D transverse Ising model, obtained by the 
2-block method. Nonzero coefficients A^r, 81,82) up to order k = 6 are listed. 
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TABLE II. Series coefficients for the dispersion E(k) /J±= J^ fe a kt „x k cos(nfc) for the singlet bound state of the Heisenberg 
ladder. Nonzero coefficients cik, n up to order k — 7 are listed. Note that the series are valid only for k > k c . 
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TABLE III. Series coefficients for dimer expansions of the energy gap E/J± of singlet bound state, triplet bound state, 
quintet antibound state, and the lower and upper edge of the continuum at k = tt for the the Heisenberg ladder. Coefficients 



of (J/J±) n up to order n = 12 are listed. 
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1.250000000 


4 


-0.203125000 


-0.476562500 


0.148437500 


-0.625000000 


-0.500000000 


5 


-0.558593750 


-0.742187500 


-0.242187500 


-1.031250000 


-1.843750000 


6 


-0.356445313 


-0.399414063 


-0.198242188 


-0.595703125 


-1.119140625 


7 


0.440856934 


0.444519043 


0.219665527 


0.648925781 


1.613769531 


8 




1.282394409 


0.294692993 


1.615997314 


3.436676025 


9 




0.964994431 


-0.865842819 


1.012023926 


1.011138916 


10 




-1.139695843 


-3.052285552 


-1.200890859 


-4.719360987 


11 




-3.099767812 


-3.914894695 


-2.788565993 


-6.971628388 


12 




-1.480682586 


0.070329791 


-0.814231584 


0.478638977 



TABLE IV. The critical point (pole) and critical index (residue) obtained from [n/m] Dlog Pade approximants to the series 
for the binding energy at k — 3%/ 5 for the triplet bound state of the Heisenberg ladder. An asterisk denotes a defective 



approximant. 



n 




[(n - 2)/n] 


[(n-l)/n] 


[n/n] 


[(« + !)/"] 


[(n + 2)/n] 






pole (residue) 


pole (residue) 


pole (residue) 


pole (residue) 


pole (residue) 


ri- 


2 




0.13342(2.100484) 


0.13067(1.917138)* 


0.13163(1.993828) 


0.13173(2.002998) 


ii— 


3 


0.13169(1.999149) 


0.13172(2.001628) 


0.13170(2.000083) 


0.13170(2.000146) 


0.13170(1.999905) 


ri- 


4 


0.13171(2.000658) 


0.13170(2.000143) 


0.13170(2.000097) 


0.13170(1.999987) 


0.13170(1.999999) 


ii — 


5 


0.13170(1.999826)* 


0.13170(1.999969) 


0.13170(1.999997) 






ri— 


6 


0.13170(2.000016) 
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